<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB R2018b"><title>Stiffness matrix of rectangular QUAD finite elements</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.7999992370605px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S2 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S3 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px;  }
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px;  }
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px;  }
.S6 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234001159668px;  }
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px;  }
.S8 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234001159668px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span style=' font-weight: bold;'>Stiffness matrix of rectangular QUAD finite elements</span></h1><div  class = 'S1'><span>We want to evaluate Stiffness matrix of a rectangular a x Gamma*a QUAD element</span></div><div  class = 'S1'><span></span></div><div  class = 'S1'><img src = "" width = "450" height = "192" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S1'><span>Ke(1,1)=(1/(1-nu^2) *E*h * (3-nu)/6</span></div><div  class = 'S1'><span>--&gt; (3-nu)*4 --&gt; A11(1)=12, B11(1)=-4 etc...</span></div><div  class = 'S1'><img src = "" width = "593" height = "102" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S1'><span>need to multiply by h in top88</span></div><h2  class = 'S2'><span>Use 2D linear elasticity</span></h2><div  class = 'S1'><span> using the Stress-strain relation for isotropic material</span></div><div  class = 'S1'><img src = "" width = "457" height = "100" alt = "" style = "vertical-align: baseline"></img></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: pre;"><span>clear </span><span class="warning_squiggle_rte" style="color: rgb(160, 32, 240);">all</span><span>; close </span><span style="color: rgb(160, 32, 240);">all</span><span>; </span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre;"><span>syms </span><span style="color: rgb(160, 32, 240);">E nu real</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: pre;"><span>C_sigma =  ([E/(1-nu^2), E*nu/(1-nu^2), sym(0); E*nu/(1-nu^2), E/(1-nu^2), sym(0); sym(0), sym(0), 0.5*E*(1-nu)/(1-nu^2)])</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="81AC0170" data-testid="output_0" data-width="698" data-height="146" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">C_sigma =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -53px;"><img src="" width="213" height="118"></span></span></div></div></div></div></div><h2  class = 'S2'><span>Shape functiton for a rectangular a x gamma*a QUAD element</span></h2><div  class = 'S1'><span>Shape functions for rectangular elements are product of Lagrange interpolations in the two coordinate directions. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: pre;"><span>syms </span><span style="color: rgb(160, 32, 240);">a gamma x y real</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'></div></div><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: pre;"><span>N_1=1/4 * (1 - (2 * x)/a) * (1 - (2 * y)/(a * gamma))</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="1AEE4F03" data-testid="output_1" data-width="698" data-height="63" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">N_1 =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -15px;"><img src="" width="129.5" height="35"></span></span></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre;"><span>N_2=1/4 * (1 + (2 * x)/a) * (1 - (2 * y)/(a * gamma))</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="9AC79937" data-testid="output_2" data-width="698" data-height="63" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">N_2 =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -15px;"><img src="" width="140" height="35"></span></span></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre;"><span>N_3=1/4 * (1 + (2 * x)/a) * (1 + (2 * y)/(a * gamma))</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="58C0FF6E" data-testid="output_3" data-width="698" data-height="63" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">N_3 =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -15px;"><img src="" width="129.5" height="35"></span></span></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre;"><span>N_4=1/4 * (1 - (2 * x)/a) * (1 + (2 * y)/(a * gamma))</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="880287E1" data-testid="output_4" data-width="698" data-height="63" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">N_4 =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -15px;"><img src="" width="140" height="35"></span></span></div></div></div></div></div><div  class = 'S8'><img src = "" width = "460" height = "241" alt = "" style = "vertical-align: baseline"></img></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: pre;"><span>N=[N_1 sym(0) N_2 sym(0) N_3 sym(0) N_4 sym(0) ;</span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: pre;"><span>   sym(0) N_1 sym(0) N_2 sym(0) N_3 sym(0) N_4]</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="BABAF579" data-testid="output_5" data-width="698" data-height="355" data-hashorizontaloverflow="false" style="width: 728px; max-height: 366px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">N =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -158px;"><img src="" width="209" height="327"></span></span></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre;"><span>N_check=subs(N,[a gamma], [1 1])</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="3B6D41A9" data-testid="output_6" data-width="698" data-height="355" data-hashorizontaloverflow="false" style="width: 728px; max-height: 366px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">N_check =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -158px;"><img src="" width="209" height="327"></span></span></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre;"><span>Bx= </span><span class="warning_squiggle_rte">[</span><span>diff(N_check(1,:), x)]</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="469A68FF" data-testid="output_7" data-width="698" data-height="64" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">Bx =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -15px;"><img src="" width="274" height="36"></span></span></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre;"><span>By= </span><span class="warning_squiggle_rte">[</span><span>diff(N_check(2,:), y)]</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="781B5717" data-testid="output_8" data-width="698" data-height="64" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">By =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -15px;"><img src="" width="272.5" height="36"></span></span></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: pre;"><span>Bxy= </span><span class="warning_squiggle_rte">[</span><span>diff(N_check(2,:), x)+diff(N_check(1,:), y)]</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="19728433" data-testid="output_9" data-width="698" data-height="64" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">Bxy =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -15px;"><img src="" width="381.5" height="36"></span></span></div></div></div></div></div><div  class = 'S8'><span>The stain displacement matrix can be computed analytically using diff</span></div><div  class = 'S1'><img src = "" width = "483" height = "180" alt = "" style = "vertical-align: baseline"></img></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: pre;"><span>B=[diff(N_1, x), sym(0), diff(N_2, x), sym(0), diff(N_3, x), sym(0), diff(N_4, x), sym(0); </span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre;"><span>    sym(0), diff(N_1, y), sym(0), diff(N_2, y), sym(0), diff(N_3, y), sym(0), diff(N_4, y); </span><span style="color: rgb(0, 0, 255);">...</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: pre;"><span>    diff(N_1, y), diff(N_1, x), diff(N_2, y), diff(N_2, x), diff(N_3, y), diff(N_3, x), diff(N_4, y), diff(N_4, x)]</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="B75A9384" data-testid="output_10" data-width="698" data-height="451" data-hashorizontaloverflow="false" style="width: 728px; max-height: 462px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">B =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -206px;"><img src="" width="254" height="423"></span></span></div></div></div></div></div><div  class = 'S8'><span>The stiffness matrix integrand is then (h is the thickness):</span></div><div  class = 'S1'><img src = "" width = "540" height = "108" alt = "" style = "vertical-align: baseline"></img></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: pre;"><span>kin=B'*C_sigma*B; </span><span style="color: rgb(34, 139, 34);">%to integrate twice</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: pre;"><span>KE = int(int(kin, x, -(a/2), a/2), y, -((a * gamma)/2), (a * gamma)/2);</span></span></div></div></div><h2  class = 'S2'><span>Substitution with 1x1 quad (for topology optimization)</span></h2><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: pre;"><span>KE=subs(KE, [a gamma], [1 1])</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement scrollableOutput" uid="2C3BE727" data-testid="output_11" data-width="698" data-height="811" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">KE =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -386px;"><img src="" width="315" height="783"></span></span></div></div></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre;"><span style="color: rgb(34, 139, 34);">%check KE(1,1)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: pre;"><span>KE(1,1)</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsSymbolicElement" uid="366EB20D" data-testid="output_12" data-width="698" data-height="66" data-hashorizontaloverflow="false" style="width: 728px; max-height: 261px;"><div class="symbolicElement"><div class="embeddedOutputsVariableElement">ans =&nbsp;</div><span class="MathEquation displaySymbolicElement" style="font-size: 15px;"><span style="vertical-align: -16px;"><img src="" width="64" height="38"></span></span></div></div></div></div></div><div  class = 'S8'><span>Right ? Ke(1,1)=((3-nu)/6*(1-nu^2)) *E --&gt; for SI UNITS please multiply by h</span></div></div><br>
<!-- 
##### SOURCE BEGIN #####
%% *Stiffness matrix of rectangular QUAD finite elements*
% We want to evaluate Stiffness matrix of a rectangular a x Gamma*a QUAD element
% 
% 
% 
% 
% 
% Ke(1,1)=(1/(1-nu^2) *E*h * (3-nu)/6
% 
% REPLACE_WITH_DASH_DASH> (3-nu)*4 REPLACE_WITH_DASH_DASH> A11(1)=12, B11(1)=-4 etc...
% 
% 
% 
% need to multiply by h in top88
%% Use 2D linear elasticity
%  using the Stress-strain relation for isotropic material
% 
% 

clear all; close all; 
syms E nu real
C_sigma =  ([E/(1-nu^2), E*nu/(1-nu^2), sym(0); E*nu/(1-nu^2), E/(1-nu^2), sym(0); sym(0), sym(0), 0.5*E*(1-nu)/(1-nu^2)])
%% Shape functiton for a rectangular a x gamma*a QUAD element
% Shape functions for rectangular elements are product of Lagrange interpolations 
% in the two coordinate directions. 

syms a gamma x y real

N_1=1/4 * (1 - (2 * x)/a) * (1 - (2 * y)/(a * gamma))
N_2=1/4 * (1 + (2 * x)/a) * (1 - (2 * y)/(a * gamma))
N_3=1/4 * (1 + (2 * x)/a) * (1 + (2 * y)/(a * gamma))
N_4=1/4 * (1 - (2 * x)/a) * (1 + (2 * y)/(a * gamma))
%% 
% 

N=[N_1 sym(0) N_2 sym(0) N_3 sym(0) N_4 sym(0) ;...
   sym(0) N_1 sym(0) N_2 sym(0) N_3 sym(0) N_4]
N_check=subs(N,[a gamma], [1 1])
Bx= [diff(N_check(1,:), x)]
By= [diff(N_check(2,:), y)]
Bxy= [diff(N_check(2,:), x)+diff(N_check(1,:), y)]
%% 
% The stain displacement matrix can be computed analytically using diff
% 
% 

B=[diff(N_1, x), sym(0), diff(N_2, x), sym(0), diff(N_3, x), sym(0), diff(N_4, x), sym(0); ...
    sym(0), diff(N_1, y), sym(0), diff(N_2, y), sym(0), diff(N_3, y), sym(0), diff(N_4, y); ...
    diff(N_1, y), diff(N_1, x), diff(N_2, y), diff(N_2, x), diff(N_3, y), diff(N_3, x), diff(N_4, y), diff(N_4, x)]
%% 
% The stiffness matrix integrand is then (h is the thickness):
% 
% 

kin=B'*C_sigma*B; %to integrate twice
KE = int(int(kin, x, -(a/2), a/2), y, -((a * gamma)/2), (a * gamma)/2);
%% Substitution with 1x1 quad (for topology optimization)

KE=subs(KE, [a gamma], [1 1])
%check KE(1,1)
KE(1,1)
%% 
% Right ? Ke(1,1)=((3-nu)/6*(1-nu^2)) *E REPLACE_WITH_DASH_DASH> for SI UNITS please multiply 
% by h
##### SOURCE END #####
--></body></html>